The performance of the quantum adiabatic algorithm 
on random instances of two optimization problems on regular hypergraphs 
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In this paper we study the performance of the quantum adiabatic algorithm on random instances 
of two combinatorial optimization problems, 3-regular 3-XORSAT and 3-regular Max-Cut. The 
cost functions associated with these two clause-based optimization problems are similar as they are 
both defined on 3-regular hypergraphs. For 3-regular 3-XORSAT the clauses contain three variables 
and for 3-regular Max-Cut the clauses contain two variables. The quantum adiabatic algorithms 
we study for these two problems use interpolating Hamiltonians which are stoquastic and therefore 
amenable to sign-problem free quantum Monte Carlo and quantum cavity methods. Using these 
techniques we find that the quantum adiabatic algorithm fails to solve either of these problems 
efficiently, although for different reasons. 
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I. INTRODUCTION 

The Quantum Adiabatic Algorithm (QAA) [l[ is an al- 
gorithm for solving optimization problems using a quan- 
tum computer. The optimization problem to be solved 
is defined by a cost function which acts on N bit strings. 
The computational task is to find the global minimum of 
the cost function. 

To use the QAA, the cost function is first encoded in a 
quantum Hamiltonian Hp (called the 'problem Hamilto- 
nian') that acts on the Hilbert space of N spin | particles. 
The problem Hamiltonian is written as a function of a z 
Pauli-matrices and is therefore diagonal in the computa- 
tional basis. The ground state of Hp corresponds to the 
solution (i.e., lowest cost bit string) of the optimization 
problem. 

To find the ground state of the problem Hamiltonian, 
the system is first prepared in the ground state of an- 
other Hamiltonian Hg, known as the beginning Hamil- 
tonian. The beginning Hamiltonian does not commute 
with the problem Hamiltonian and must be chosen so 
that its ground state is easy to prepare. Here we use the 
standard choice 



tween the two Hamiltonians 



H(s) = (1 - s)H B + sH P . 



(1) 



where s(t) is a parameter varying smoothly with time, 
from s(0) = to s(T) = 1 at the end of the algorithm 
after a total evolution time T ■ 

If the parameter s(t) is changed slowly enough, the adi- 
abatic theorem of Quantum Mechanics (sec. e.g., Rcfs. [2] 
and (H) ensures that the system will stay close to the 
ground state of the instantaneous Hamiltonian through- 
out the evolution. After time T the state obtained will be 
close to the ground state of Hp. A final measurement of 
the state in the Pauli z basis then produces the solution 
of the optimization problem. 

The runtime T must be chosen to be large enough so 
that the adiabatic approximation holds: this condition 
determines the efficiency, or complexity, of the QAA. A 
condition on T can be given in terms of the eigenstates 
{|m)} and eigenvalues {E m } of the Hamiltonian H(s), 
as HH 



|max s Vi (s)| 
(A£ min ) 2 



(2) 



H B = 



N 

E 



which has a product state as its ground state. 

The Hamiltonian of the system is slowly modified from 
Hb to Hp. Here we consider a linear interpolation be- 



where A£" m i n is the minimum of the first excitation gap 
A-E m i n = min s Ai? with AE = E\ — Eq, and V m Q = 
(0\dH/ds\m). 

Typically, matrix elements of H(s) scale as a low poly- 
nomial of the system size A, and the question of whether 
the runtime is polynomial or exponential as a function of 
N therefore depends on how the minimum gap A_E m j n 
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scales with N. If the gap becomes exponentially small 
at any point in the evolution, then the computation re- 
quires an exponential amount of time and the QAA is 
inefficient. The dependence of the minimum gap on the 
system size for a given problem is therefore a central issue 
in determining the complexity of the QAA. 

A notable feature of the interpolating Hamiltonian (fT]) 
is that it is real and all of its off diagonal matrix elements 
are non-positive. Hamiltonians which have this property 
are called stoquastic [f|. There is complexity-theoretic 
evidence that some computational problems regarding 
the ground states of stoquastic Hamiltonians are easier 
than the corresponding problems for more general Hamil- 
tonians Q. It may be the case that quantum adiabatic 
algorithms using stoquastic interpolating Hamiltonians 
(such as the ones we consider here) are no more powerful 
than classical algorithms-this remains an intriguing open 
question. 

An interesting question about the QAA is how it per- 
forms on "hard" sets of problems - those for which all 
known algorithms take an exponential amount of time. 
While early studies of the QAA done on small systems 
(N < 24) Q indicated that the time required to solve 
one such problem might scale polynomially with N, sev- 
eral later studies using larger system sizes gave evidence 
that this may not be the case. 

References pi Q show that adiabatic algorithms will 
fail if the initial Hamiltonian is chosen poorly. Recent 
work has elucidated a more subtle way in which the adi- 
abatic algorithm can fail [Tol - |l^ | . The idea of these works 
is that a very small gap can appear in the spectrum of 
the interpolating Hamiltonian due to an avoided crossing 
between the ground state and another level correspond- 
ing to a local minimum of the optimization problem. The 
location of these avoided crossings moves towards s = 1 
as the system size grows. They have been called "per- 
turbative crosses" because it is possible to locate them 
using low order perturbation theory. Altshuler et al. [l3| 
have argued that this failure mode dooms the QAA for 
random instances of NP-complete problems. However, 
the arguments of Altshuler et al. have been criticized by 
Knysh and Smelyanskiy (TH |. 

Young et al. [l6l , [ItJ recently examined the perfor- 
mance of the QAA on random instances of the constraint 
satisfaction problem called l-in-3 SAT (to be described 
in the next section) and showed the presence of avoided 
crossings associated with very small gaps. These 'bot- 
tlenecks' appears in a larger and larger fraction of the 
instances as the problem size N increases, indicating the 
existence of a first order quantum phase transition. This 
leads to an exponentially small gap for a typical instance, 
and therefore also to the failure of adiabatic quantum op- 
timization. 

It is not yet clear to what extent the above behavior 
found for l-in-3 SAT is general and whether it is a fea- 
ture inherent to the QAA that will plague most if not 
all problems fed into the algorithm or something more 
benign than this. Previous work (T^-HU had argued that 



a first order quantum phase transition occurs for a broad 
class of random optimization models. 

In this paper we contrast the performance of the quan- 
tum adiabatic algorithm on random instances of two com- 
binatorial optimization problems. The first problem we 
consider is 3-XORSAT on a random 3-regular hyper- 
graph, which was studied previously in Ref. [jgj . In- 
terestingly, although this computational problem is clas- 
sically easy-an instance can always be solved in poly- 
nomial time on a classical computer by using Gaussian 
elimination— it is known that classical algorithms that do 
not use linear alg ebra are stymied by this problem (22| - 
Hi|. In Ref. pjf it was shown that the QAA fails to 
solve this problem in polynomial time. In this paper we 
provide more numerical evidence for this. We also fur- 
nish a duality transformation that helps to understand 
properties of this model. 

The second computational problem we consider is 
Max-Cut on a 3-regular graph. This problem is NP- 
hard. However we consider random instances, for which 
the computational complexity is less well understood. 

A nice feature of these problems is that the regularity 
of the associated hypergraphs constrains the two ensem- 
bles of random instances. Studying the performance of 
the QAA for these problems, we therefore expect to see 
smaller instance-to-instance differences than for the un- 
constrained ensembles of instances. 

We use two different methods to study the performance 
of the QAA. The first method is quantum Monte Carlo 
simulation. It is a numerical method that is based on 
sampling paths from the Taylor expansion of the parti- 
tion function of the system. Using this method we can 
extract, for a given instance, the thermodynamic proper- 
ties (in particular the ground state energy) as well as the 
eigenvalue gap for the interpolating Hamiltonian H(s). 
This allows us to investigate the size dependence of the 
typical minimum gap of the problem from which we can 
extrapolate the large-size scaling of the computation time 
T of the QAA. The second approach is a quantum cav- 
ity method. It is a semi-analytical method that allows 
us to compute the thermodynamic properties averaged 
over the ensemble of instances in the limit N — > oo. It 
leads to a set of self-consistent equations that can be 
solved analytically in some classical examples (26l . [27j . 
However in the quantum case the equations are more 
complicated and are solved numerically [2^, (29|. The 
method is not exact on general graphs. For locally tree- 
like random graphs, it provides the exact solution of the 
problem if some assumptions on the Gibbs measure are 
satisfied (2(| [2?], H(j ■ As we will discuss below the cavity 
method we use in this paper gives the exact result for 
3-XORSAT, while it only gives an approximation for the 
Max-Cut problem. 

Using these methods we conclude that the quantum 
adiabatic algorithm fails to solve both problems effi- 
ciently, although in a qualitatively different way. 

The plan of this paper is as follows. In Section [XT] 
we describe the two computational problems that we in- 
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vestigate. In Sec. IIIII we discuss the methods that we 
use to obtain our results. These results are presented in 
Sec. IIVI and our conclusions are summarized in Sec. [V] 
Some parts of this paper have previously appeared in the 
PhD thesis of one of the authors [1H . 



II. MODELS 

We now discuss in detail the two computational prob- 
lems 3- regular 3-XORSAT and 3- regular Max-Cut. 

When studying the efficiency of the QAA numeri- 
cally 0, [Tsl . it is convenient to consider instances 
with a unique satisfying assignment (USA) for reasons 
that will be explained in Sec. IIII Al On the other hand, 
the quantum cavity method is designed to study the en- 
semble of random instances with no restrictions on the 
number of satisfying assignments. In this section we spec- 
ify the random ensembles of instances that we investigate 
in this paper. 



A. 3-regular 3-XORSAT 

The 3-XORSAT problem is a clause based constraint 
satisfaction problem. An instance of such a constraint 
satisfaction problem is specified as a list of M logical 
conditions (clauses) on a set of N binary variables. The 
problem is to determine whether there is an assignment 
to N bits which satisfies all M clauses. 

In the 3-XORSAT problem each clause involves three 
bits. A given clause is satisfied if the sum of the three bits 
(mod 2) is a specified value (either or 1, depending on 
the clause). We consider the "3-regular" case where every 
bit is in exactly three clauses which implies M = N. This 
model has already been considered by Jorg et al. [l9j . 
The factor graph for an instance of 3-regular 3-XORSAT 
is sketched in Fig. [TJ 

Since this problem just involves linear constraints (mod 
2), the satisfiability problem can be solved in polyno- 
mial time using Gaussian elimination. However, it is well 
known that this problem presents difficulties for solvers 
that do not use linear algebra (see, e.g. Refs. [22T - [25l |). 

We associate each instance of 3-regular 3-XORSAT 
with a problem Hamiltonian Hp that acts on N spins. 
Each clause is mapped to an operator which acts nontriv- 
ially on the spins involved in the clause. The operator 
for a given clause has energy zero if the clause is satisfied 
and energy equal to I if it is not, so 

H P = J2 ( I 3 ° a *°> °* j . (3) 

Here each clause c € {1,...,A} is associated with the 3 
bits *i c)*2o*3c an d a coupling J c € {±1} which tells us 
if the sum of the bits mod 2 should be or 1 when the 
clause is satisfied. 




Figure 1: Factor graph of a small part of an instance of the 
3-regular 3-XORSAT problem. In the full factor graph, each 
clause (□) is connected to exactly three bits (O) an d each bit 
is connected to exactly three clauses, so there are no leaves 
and the graph closes up on itself. 

1. Random Instances of 3-regular 3-XORSAT 

As inRef. [H], we consider both the random ensemble 
of instances of this problem and the random ensemble 
of instances which have a unique satisfying assignment 
(USA). In the 3-XORSAT problem as N -> oo, instances 
with a USA are a nonzero fraction, about 0.285 [l9| . of 
the set of all instances, so the random ensemble of USA 
instances should be a good representation of the fully 
random ensemble. 

All satisfiable instances (and in particular instances 
with a USA) have the property that the cost function 
Eq. (J3j) can be mapped unitarily into the form 

. , / 1 _ _tl,C *3, c _.i 3 ,c\ 

#p = E f Z 2 Z °' j (4) 
by a product of bit flip operators. 



2. Previous Work 

Reference [l9| studied the performance of the QAA 
on the random ensemble of instances of 3-regular 3- 
XORSAT using quantum cavity method and quantum 
Monte Carlo simulation. They also studied the ensemble 
of random instances with a USA using exact numerical 
diagonalization. This work gave evidence that there is 
a first order quantum phase transition which occurs at 
s c ~ | in the ground state. Their results also demon- 
strate that the minimum gap is exponentially small as a 
function of N at the transition point. 
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3. Duality Transformation 



Note that 



In this section we demonstrate a duality mapping 
for the ensemble of random instances of 3-regular 3- 
XORSAT with a unique satisfying assignment. This du- 
ality mapping explains the critical value = | of the 
quantum phase transition in this model [19|. Consider 
the Hamiltonian 



JV 



1 - <r\ 



N 



1 - al 



(5) 

Here, the first term is the beginning Hamiltonian and the 
second term is the problem Hamiltonian for an instance 
of 3-regular 3-XORSAT with a unique satisfying assign- 
ment. The 3-regular hypergraph specifying the instance 
can be represented by a matrix M where 



Ma = 



1 , if bit j is in clause i 
, otherwise. 



and where M has 3 ones in each row and 3 ones in each 
column. The fact that there is a unique satisfying as- 
signment 000... is equivalent to the statement that the 
matrix M is invertible over F^. To see this, consider the 
equation (with addition mod 2) 



Mv = 0. 

This equation has the unique solution v — if and only 
if there is a unique satisfying assignment for the given 
instance. This is also the criterion for the matrix M to 
be invertible. 

The duality that we construct shows that the spectrum 
of -ff(s) is the same as the spectrum of i?DUAL(l — s) 
where Houal is obtained by replacing the problem 
Hamiltonian hypergraph by its dual-that is to say, the in- 
stance corresponding to a matrix M is mapped to the in- 
stance associated with M T . The ground state energy per 
spin (averaged over all 3-regular instances with a unique 
satisfying assignment) is symmetric about s = ^ and the 
first order phase transition observed in Ref. |19| occurs 
at s = h. For each c = 1, N define the operator 



X c = a z " cr/' . (6) 
We also define, for each clause c, a bit string if 
F = M~ x e c . 

Here e c is the unit vector with components (e c )i = Si C . 
Note that y° is the unique bit string which violates clause 
c and satisfies all other clauses. Such a bit string is guar- 
anteed to exist since M is invertible. Let yf denote the 
ith bit of the string y^. Define, for each c = 1, N, 



N 



z=R [4] 



(7) 



and 



{Z c ,X c }=0 



[Z C ,X C ,] = for c^c'. 



For each bit i = 1,...,N let ci(i), 02(1), 03(1) be the 
clauses which bit i participates in. Then 



(8) 



This follows from the fact that 



and so 



Met = e Cl (j) + e C2 (i) + e C3 ^ 

&i = (e Cl (i) +e C2(l) +e C3(l) ) 

= y 01 ^ + y 02 ^ + y° 3 ^ . 

The above equation and the definition Eq. (0 show 
Eq. ©. Now using Eqs. © and © write 



N 



H(s) = (l-s)Y^ 



1 - o\ 



N 



1 - o\^ c a?' c o\^ c 



N 



1 - Z c 1 {i) Z c 2 {i) Z c 3 {i) 



N 



l-Xr 



(9) 
(10) 

(11) 
(12) 



The X and Z operators satisfy the same commutation 
relations as the operators a x and a z . Comparing Eq. Q12p 
with Eq. §5§ we conclude that the spectrum of H(s) is 
the same as the spectrum of Hdual{1 — s). 

In Fig. [5] we show the first four energy levels of the 
interpolating Hamiltonian H(s) = sHb + (1 — s)Hp as a 
function of s for one 16-bit instance of 3-XORSAT. The 
duality transformation means that these energy levels are 
the same as for the interpolating Hamiltonian H'(s) = 
sHp^dual + (1 — s)Hb which involves the dual instance. 
Evident from the figure is the apparent symmetry of the 
energy levels around s = 1/2. In this case the instance 
and its dual are similar from the point of view of the 
QAA. 

The duality argument given here has implications for 
the phase transition which occurs in the ensemble of ran- 
dom instances of 3-regular 3-XORSAT as N — s- 00. Our 
numerics in section HVl show that for large N , the ground 
state energy per spin as a function of s (averaged over 
the ensemble of instances) has a nonzero derivative as 
s — > i. The duality transformation given here implies 
that this curve is symmetric about s = 5. So there is 
a discontinuity in the derivative of this curve at s = 5, 
which is associated with a first order phase transition. 
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Figure 2: (Color online) First four energy levels of the in- 
terpolating Hamiltonian for a 16-bit instance of the 3-regular 
3-XORSAT problem. The energy curves for this instance are 
close to being symmetric about s — 1/2. Our duality transfor- 
mation means that sending s —¥ (1— s) we obtain the spectrum 
of the interpolating Hamiltonian for a different instance from 
the same ensemble, obtained by interchanging the clauses and 
bits. 

B. 3-regular Max-Cut 

The second model we discuss is also a clause based 
problem. The instances we consider are not satisfiable 
and we are interested in finding the assignment which 
gives the maximum number of satisfied clauses. We view 
this problem as minimizing a cost function that computes 
the number of unsatisfied clauses. The 3-regular Max- 
Cut problem is defined on N bits, and each bit appears in 
exactly three clauses. Each clause involves two bits and 
is satisfied if and only if the sum of the two bits (modulo 
2) is 1. The number of clauses is therefore M = 3JV/2. 
The problem Hamiltonian is 

Hp ^(l±ivpvfy (13) 

The ground state of this Hamiltonian encodes the solu- 
tion to the Max-Cut problem. 

The model can also be viewed as an antiferromagnet 
on a 3-regular random graph. Because the random graph 
in general has loops of odd length, it is not possible to 
satisfy all of the clauses. 

The Max-Cut problem is NP-hard and accordingly 
there is no known classical polynomial time algorithm 
which computes the ground state energy of the problem 
Hamiltonian (|13[) . Indeed, even achieving a certain ap- 
proximation to the ground state energy is hard, which 
follows from the fact that it is NP hard to approximate 
the Max-Cut of 3-regular graphs to within a multiplica- 



tive factor 0.997 [Hj]. Interestingly, however, there is a 
classical polynomial time algorithm which achieves an 
approximation ratio of at least 0.9326 [33[. 

1. Random Instances of 3-regular Max-Cut 

Using the quantum cavity method we study the en- 
semble of random instances of 3-regular Max-Cut. 

The random instances we studied using quantum 
Monte Carlo simulation were restricted to those which 
have exactly 2 minimal energy states (note that this is the 
smallest number possible since the problem is symmetric 
under flipping all the spins) and for which the ground 
state energy of the problem Hamiltonian is equal to \N . 
We choose to study instances with a unique satisfying 
assignment (up to the bit-flip symmetry of this problem) 
because it is numerically more convenient for the extrac- 
tion of the relevant gap (to the first even state) . For the 
range of sizes studied, was found numerically to be 
the most probable value of the ground state energy. The 
restriction to instances with a fixed Max-Cut (|-/V~) fur- 
ther reduces the instance-to-instance fluctuations. How- 
ever, this choice affects the ensemble averaged value of 
thermodynamic observables (e.g. the average energy of 
fully random instances is different from N/8), making it 
more difficult to compare the quantum Monte Carlo re- 
sults with our quantum-cavity results on the fully random 
ensemble. We expect (and find numerically) that this set 
of instances makes up an exponentially small fraction of 
the whole random ensemble for large N. 

2. Previous work 

Laumann et al. [28[ used the quantum cavity method 
to study the transverse field spin glass with the problem 
Hamiltonian 

^=E( 1+ T a? " ) (14) 

where each J c is chosen to be +1 or —1 with equal prob- 
ability. 

In general there is no "gauge transformation" equiv- 
alence between this problem Hamiltonian and the anti- 
ferromagnetic problem Hamiltonian Eq. (|13p . However 
we do expect these models to exhibit similar properties 
since a random graph is locally tree-like, and on a tree 
such a gauge transformation does exist, see Ref. [34| for 
a discussion of this point in the case where there is no 
transverse field present. 

Laumann et al. found that this system exhibits a sec- 
ond order phase transition as a function of the transverse 
field. Their method is similar to the quantum cavity 
method that we use, although the numerics performed 
in Ref. [28| have some systematic errors which our cal- 
culations avoid. The method used in Ref. [28[ is a dis- 
crete imaginary time formulation of the quantum cavity 



6 



method which has nonzero Trotter error, whereas our cal- 
culation works in continuous imaginary time [29| where 
this source of error is absent. Our calculation also does 
not use the approximation used in Ref. [28[ where the "ef- 
fective action" of a path in imaginary time is truncated 
at second order in a cluster expansion. 

III. METHOD 
A. Quantum Monte Carlo 

The complexity of the QAA algorithm is determined by 
the size dependence of the "typical" minimum gap of the 
problem. Following Refs. (l^, [l7l . l2ll |. we analyze the size- 
dependence of these gaps by considering (typically) 50 
instances for each size, and then extracting the minimum 
gap for each of them. For each instance, we perform 
quantum Monte Carlo simulations for a range of s values 
and hunt for the minimum gap. We then take the median 
value of the minimum gap among the different instances 
for a given size to obtain the "typical" minimum gap. 

Quantum Monte Carlo simulation works by sampling 
random variables from a probability distribution (over 
some configuration space) which contains information 
about the quantum system of interest. The probability 
distribution is sampled by Markov chain Monte Carlo, 
and properties of the quantum system to be studied 
are obtained as expectation values. Different quantum 
Monte Carlo methods are based on different ways of as- 
sociating probability distributions to a quantum system. 

In our simulations we use a quantum Monte Carlo tech- 
nique known as the stochastic series expansion (SSE) al- 
gorithm [35|, [3(|. In this method, the probability dis- 
tribution associated with the quantum system is derived 
from the Taylor series expansion of the partition function 
Tr[e -,3ff ] at inverse temperature /?. This is in contrast to 
other quantum Monte Carlo techniques which are based 
on the path integral expansion of the partition function. 
Whereas some of these techniques have systematic errors 
because the path integral expansions used are inexact, 
the SSE that we use has no such systematic error. 

A second feature of the SSE method that we use is 
that the Markov chain used to sample configurations al- 
lows global, as well as local, updates, which leads to faster 
equilibration. We further speed up equilibration by im- 
plementing "parallel tempering" [37] . where simulations 
for different values of s are run in parallel and spin con- 
figurations with adjacent values of s are swapped with 
a probability satisfying the detailed balance condition. 
Traditionally, parallel tempering is performed for systems 
at different temperatures, but here the parameter s plays 
the role of (inverse) temperature (38| . 

The details of our implementations of the SSE method 
for 3-XORSAT and Max-Cut are slightly different and 
are further discussed in Appendix[3] Moreover, the Max- 
Cut problem can not be simulated with the 'traditional' 
SSE method because of its symmetry under flipping all of 



the spins: For a given instance of Max-Cut, every eigen- 
state of the interpolating Hamiltonian H (s) is either even 
or odd under this symmetry transformation. Since the 
ground state is even, here we are interested in the eigen- 
value gap to the first even excited state. We therefore 
design our quantum Monte Carlo simulation so that it 
works in the subspace of even states. The modified algo- 
rithm is detailed in Appendix [Bj 

In our simulations we extract the gap from imaginary 
time-dependent correlation functions. The gap of the 
system for a given instance and a given s value is ex- 
tracted by analyzing measurements of (imaginary) time- 
dependent correlation functions of the type 

C A (r) = (i(r)i(O)) - (Af , (15) 

where the operator A is some measurable physical quan- 
tity. It is useful to optimize the choice of correlation 
functions such that the contribution from the first excited 
state, to = 1 in Eq. (|16p below, is as large as possible rela- 
tive to the contributions from higher excited states. One 
way of doing this, which was used in some of the runs, is 
described in Ref. [39j . 

The evaluation of (A) 2 in the above equation is com- 
puted from the product (A)^ (A)^ where the two in- 
dices correspond to different independent simulations of 
the same system. This eliminates the bias stemming from 
straightforward squaring of the expectation value. 

In the low temperature limit, T AEi where AE\ = 
Ei — Eq 1 the system is in its ground state so the 
imaginary-time correlation function is given by 

C A {r) = £ |(0|i|m)| 2 (e-^ + e -A^-r)) , 

m— 1 

(16) 

where AE m = E m — Eo. At long times, r, the correlation 
function is dominated by the smallest gap AEi (as long 
as the matrix element |(0|^4|1)| 2 is nonzero). On a log- 
linear plot Ca(t), then has a region where it is a straight 
line whose slope is the negative of the gap. This can 
therefore be extracted by linear fitting. A more detailed 
description of the method may be found in Ref. [39j . 

B. The quantum cavity method 

The quantum cavity method [28|,[29| is a technique that 
is used to study thermodynamic properties of transverse 
field spin Hamiltonians. In our implementation we use 
the continuous imaginary time method from Ref. [29j . 
Quantum Cavity methods have now been used to study 
a number of problems including the ferromagnet on the 
Bethe lattice in uniform [29| and random [4Cj trans- 
verse field, the spin glass on the Bethe Lattice (28|, 3- 
regular 3-XORSAT [lg, and the quantum Biroli-Mezard 
model [il|. 

If the Hamiltonian is a two-local transverse field Hamil- 
tonian on a finite number of spins and if the interaction 
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graph consists of a tree (i.e if there are no loops) then 
the quantum cavity equations are exact. In this case the 
quantum cavity equations are a closed set of equations 
that exactly characterize the thermodynamic properties 
of the system at a fixed inverse temperature j3. If instead 
the interaction graph is a random regular graph with a fi- 
nite number of spins then it must have loops. As N — > oo 
we can think of it as an "infinite tree" since the typical 
size of loops in such a graph diverges. 

Quantum cavity methods (such as the one we use) for 
problems defined on random regular graphs make use of 
two properties of the system: (i) the fact that a random 
regular graph is locally tree- like; (ii) the fact that spin- 
spin correlations decay quickly as a function of distance. 
While the first property is true with probability 1 for ran- 
dom regular graphs when N — > oo, the second property 
is not always true and we now discuss it more carefully. 

The simplest case is when the Gibbs measure is char- 
acterized by a single pure state that has the clustering 
property, as in a paramagnetic phase. This happens at 
high enough temperature or large enough transverse field. 
In this case, correlations decay exponentially and the sim- 
plest version of the cavity method (the so-called "replica 
symmetric (RS)" cavity method) gives the exact result. 
Upon lowering the temperature or the transverse field, 
a phase transition towards a more complicated phase 
can be encountered. If this phase is a standard broken- 
symmetry phase (e.g. a ferromagnetic phase), then cor- 
relation decay holds provided one adds an infinitesimal 
symmetry-breaking field, and the RS cavity method still 
provides the exact result (29| . 

However, if the transition is to a spin glass phase, 
then the Gibbs measure is split into a large number of 
states, and the decorrelation property that is required by 
the cavity method only holds within each state. In this 
case, there is no explicit symmetry breaking, therefore 
the states cannot be selected by adding an infinitesimal 
external field. It turns out that refined versions of the 
cavity method must be used, that are based on assump- 
tions on the structure of these states j30|. The simplest 
assumption is that states are distributed in a uniform 
way in the phase space of the system, and leads to the 
so-called "1 step replica symmetry breaking (1RSB)" ap- 
proximation. In more complicated cases, states might 
be arranged in "clusters" leading to a hierarchical orga- 
nization, and this requires further steps of RSB (42|. A 
consistency check can be performed within the method to 
check whether a given RSB scheme gives the exact result 
or whether further RSB steps are required. 

For the XORSAT problem on random regular graphs, 
it can be rigorously shown that the 1RSB scheme gives 
the exact result in the classical case [4^, [44| , and it has 
been conjectured that the same is true for the quantum 
problem in transverse field [lj|- For the specific case 
of 3-regular 3-XORSAT investigated here, a RS calcula- 
tion is enough to get the thermodynamic properties [l9| , 
which is why we use the RS method in our simulations of 
this model. We study this problem at a temperature low 



enough that no residual temperature dependence of the 
energy is observed. Furthermore, we set the parameters 
of our calculation to be more computationally demanding 
than those of reference pj|, which allows us to achieve 
better precision. 

The study of 3-regular Max-Cut is more involved. To 
understand how well the cavity method works on this 
problem, we can look at results obtained for the classical 
3 regular spin-glass with Hamiltonian (fT4")l . We expect 
that these problems have very similar (possibly identi- 
cal) thermodynamic properties (34|. For the classical 3- 
regular spin glass it can be shown that neither the RS nor 
the 1RSB cavity method give exact results (2(1 l27l. |34| . 
and it is widely believed that an infinite number of RSB 
steps is required (this has been shown rigorously for the 
^-regular spin glass as z — > oo, which corresponds to the 
Sherrington-Kirkpatrick model [IH). However the 1RSB 
calculation gives a good approximation to the classical 
ground state energy (27} • 

To study 3-regular Max-Cut we used the 1RSB quan- 
tum cavity method with Parisi parameter m = 0. This 
level of approximation is more accurate than the replica 
symmetric approximation but less accurate than the full 
1RSB calculation. With this method we studied the 
N — > oo limit of the random ensemble of 3-regular Max- 
Cut Hamiltonians 

N 

c i—1 

We ran our calculations at finite inverse tempera- 
ture f3 = 4 and various values of the transverse field 
A. Thermodynamic expectation values with respect to 
H(X) at inverse temperature f3 can be related to ther- 
modynamic expectation values with respect to H(s) = 
(1 — s)Hb + sHp at s — and inverse temperature 
ft = M. Note that there is an s dependence introduced 
in the temperature when relating these two thermody- 
namic ensembles. 



IV. RESULTS 

In what follows we present the results of the QMC 
simulation alongside those of the quantum cavity results 
for the two problems we study here. We show that the 
QAA fails with the choice of interpolating Hamiltonians 
discussed previously; for both problems the running time 
appears to be exponentially long as a function of the 
problem size. However, the reasons for this failure are 
different for each of the models. 



A. Random 3-regular 3-XORSAT 

The 3-regular 3-XORSAT problem was studied by Jorg 
et al. who determined the minimum gap for sizes 
up to N = 24. Here, we extend the range of sizes up 
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Figure 3: (Color online) Median minimum gap as a function 
of problem size of the 3-regular 3-XORSAT problem on a 
log-linear scale. The straight-line fit is good, indicating an 
exponential dependence which in turn leads to an exponential 
complexity of the QAA for this problem. Triangles indicate 
exact-diagonalization results while the circles are the results 
of QMC simulations. 

to N = 40 by quantum Monte Carlo simulations. The 
two sets of results agree and provide compelling evidence 
for an exponential minimum gap. The duality argument 
in Sec. Ill A 31 shows that the quantum phase transition 
occurs exactly at s = s c = 1/2. Our numerics show that 
the phase transition is strongly first order, in agreement 
with Ref. [H . 

We show results for the median minimum gap as a 
function of size for the 3-regular 3-XORSAT problem in 
Fig. [3] (log-lin). A straight line fit works well for the 
log-lin plot, which provides evidence that the minimum 
gap is exponentially small in the system size. The results 
shown here generalize and agree with those obtained by 
Jorg et al. |19j . 

We also computed some ground-state properties of the 
model: the energy (H), the magnetization along the x- 
axis M x — jj X^i^i )' an d the spin-glass order parameter 
defined by: 



(17) 



These quantities, averaged over 50 instances for each size, 
are plotted in Figs. [5J [5J and [5] 

Figure [4] shows that for large system sizes differences 
between the QMC results for the ground state energy 
and the (replica symmetric) cavity results are small. Ref- 
erence [19fl has argued that the replica symmetric (RS) 
cavity method is actually exact for the thermodynamic 
properties of the 3 regular 3-XORSAT problem. To check 
this, in Fig. [7] we have expanded the vertical scale and 
show an extrapolation of the QMC results to N — oo at 
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Figure 4: (Color online) Mean energy (averaged over 50 sam- 
ple instances per size) of the 3-regular 3-XORSAT problem 
as a function of the adiabatic parameter s for different sizes 
(QMC results) compared with the RS quantum cavity calcu- 
lations. Because of the duality of the model, the true curve 
(averaged over all instances at a given value of N) is sym- 
metric about s = 1/2. The main panel shows a blowup near 
the symmetry point s — 1/2. In the inset, the entire range is 
shown. 




0.504 



Figure 5: (Color online) Magnetization along the a;-axis, 
M x = iV^ 1 ^L(crf), as a function of the adiabatic parame- 
ter s for the 3-regular 3-XORSAT problem. Results obtained 
both by QMC and the cavity method are shown. The latter 
indicates a sharp discontinuity at s = s c — 1/2. The slope 
of the QMC results at s = 1/2 increases with increasing N, 
consistent with a discontinuity at N = oo. 
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Figure 6: (Color online) The spin-glass order parameter q as 
denned in Eq. (JXTJ) as a function of the adiabatic parameter 
s for the 3-regular 3-XORSAT problem. The rapid change 
for large sizes around s =1/2 indicates a first-order quantum 
phase transition at this value of s. 



lations at s = i. In the quantum Monte Carlo data we 
see that the slope of the magnetization increases with N 
and is therefore consistent with the cavity prediction for 
N -> oo. 



B. Random 3-regular Max-Cut 

In the Monte Carlo simulations of the Max-Cut prob- 
lem we restrict ourselves to instances for which the prob- 
lem Hamiltonian has a ground state degeneracy of two, 
and for which the ground state energy is N/8. For this 
ensemble of instances, we measured the energy, the re- 
magnetization and the spin-glass order parameter using 
quantum Monte Carlo simulations. Because of the bit- 
flip symmetry of the model, we use the following different 
definition of the spin glass order parameter: 



1/2 



1 



Q = 



N(N — 1) 



z\2 



(18) 



i*7 
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0.005 0.01 0.015 
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Figure 7: (Color online) Extrapolation of the energy values 
as given by the QMC method (solid line) for different sys- 
tem sizes at s = 1/2 as compared to the value given by 
the cavity method (which is for N = oo) for the 3-regular 
3-XORSAT problem, assuming a 1/N dependence. Extrap- 
olating the QMC results to N = oo seems to given a result 
consistent with the cavity value. 

the critical value s = s c = 1/2 (where finite-size correc- 
tions are largest). The extrapolated value appears to be 
consistent with the cavity result. 

The rapid variation of M x and q shown in figures [5] 
and [5] in the vicinity of s c = 1/2 is evidence for a first- 
order transition. Figure [5] also shows a discontinuity in 
the x-axis magnetization predicted by the cavity calcu- 




0.3 0.32 0.34 0.36 0.38 0.4 0.42 0.44 0.46 
s 



Figure 8: (Color online) The spin-glass order parameter q' ob- 
tained from Monte Carlo simulations, obtained from Eq. (|18|) . 
as a function of the adiabatic parameter s for the Max-Cut 
problem. Also shown is the value of q from the cavity calcula- 
tion, which is defined differently as discussed in the text. The 
inset shows a global view over the whole range of s, indicating 
large differences between the Monte Carlo and cavity calcula- 
tions for large s. This may be due to the different ensembles 
used in the two calculations, as discussed in the text. 
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Figure 9: (Color online) Magnetization along the x direction, 
M x = N^ 1 53 (cf), as a function of the adiabatic parameter 
s for the Max-Cut problem. 



two calculations. We also note that our cavity method 
computation is performed at nonzero temperature. 



Using quantum Monte Carlo simulations we have de- 
termined the energy gap as a function of s for s in the 
range between 0.3 (i.e. well below s c ) and 0.5 (i.e. well 
above s c ) for sizes between N = 16 and 160. For the 
smaller sizes we find a single minimum in this range, 
which lies a little above s c . However, for larger sizes, we 
see a fraction of instances in which there is a minimum 
close to s c and a second, deeper, minimum for s > s c well 
inside the spin-glass phase. A set of data which shows two 
minima is presented in Fig. llH This interesting behavior 
of the minimum gaps suggests the following interpreta- 
tion: The minima found close (just above) s c correspond 
to the order-disorder quantum phase transition. Above 
s c the system is the spin-glass phase. The minima that 
are well within the spin-glass phase may correspond to 
'accidental' or perturbative crossings in the spin glass. 



In Figs . [8l [91 and [TOl we compare the QMC results with 
those of our quantum cavity method computation. Re- 
call that the cavity method results apply to the random 
ensemble of instances. Formally the value of the spin 
glass parameter q from Eq. (fT7f is zero due to the bit 
flip symmetry of the Hamiltonian. However, the cavity 
method works in the thermodynamic limit in which this 
symmetry is spontaneously broken for s greater than the 
critical value s c . For the cavity calculations, we measure 
a spin glass order parameter q which is the magnetiza- 
tion squared for each thermodynamic "state", averaged 
over the states. This becomes non-zero for s > s c ~ 0.36 
as shown in Fig. El The Monte Carlo results are consis- 
tent with this value of s c . The inset of Fig. [5] shows a 
global view of the two spin glass order parameters that 
we measure, over the whole range of s. There we see 
differences between the (different) order parameters mea- 
sured using Monte Carlo and cavity calculations for large 
s. Note that the Monte Carlo simulations take only in- 
stances with a doubly degenerate ground state for the 
problem Hamiltonian, s = 1, so q' = 1 in this limit, 
whereas the cavity calculations are done for the random 
ensemble where the instances have much larger degener- 
acy at s = 1 and so q < 1 in this limit. 

Our numerical results for the x-component of the mag- 
netization are shown in Fig. We see no evidence of 
a discontinuity in this quantity at s c ~ 0.36 for large 
N. This is in contrast with the corresponding plot for 
3-XORSAT in figure O 

Results for the energy of the Max-Cut problem, ob- 
tained both from Monte Carlo and the cavity approach 
are shown in Fig. [TOl The two agree reasonably well 
but there are differences in the spin glass phase, ,s > ,s c , 
which may be due to the different ensembles used in the 



Double-minima occurrences become more frequent as 
the system size increases. While no double minima were 
found for sizes TV = 16, 24 and 32 (within the studied 
range of s values), for sizes N = 64,128 and 160 the 
percentage of instances that exhibit such double minima 
was found to be approximately 7%, 36% and 40%, re- 
spectively (obtained from ~ 50 instances for each size). 



We have therefore performed two analyses on the data 
for the gap. In the first analysis we determine the global 
minimum (for the range of s studied) for each instance 
and determine the median over the instances. There are 
about 50 instances for each size. This data is presented 
in Table U and is plotted in Fig. [TO] both as log-lin (main 
figure) and log-log (inset). A straight line fit works well 
for the log-lin plot (goodness of fit parameter Q = 0.57), 
provided that we omit the two smallest sizes. However, 
a straight-line fit works much less well for the log-log 
plot (Q = 2.7 x 10 -3 ), again omitting the two smallest 
sizes, because the data for the largest size lies below the 
extrapolation from smaller sizes. If smaller points do not 
lie on the fit, it is possible that the fit is correct and the 
deviations are due to corrections to scaling. However, 
if the largest size shows a clear deviation then the fit 
can not describe the asymptotic large- ./V behavior. From 
these fits we conclude that an exponentially decreasing 
gap is preferred over a polynomial gap. 
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Figure 10: (Color Online) Energy as a function of the adi- 
abatic parameter s for the Max-Cut problem. The cavity 
results computed at inverse temperature /3 = - are depicted 
by the dashed line. The lower panel is a blow up of the region 
around the maximum, which illustrates the difference between 
the Monte Carlo and cavity results. Some of this difference 
may be due to the different ensembles used, as discussed in 
the text. 



Table I: Median minimum gap for 3-regular Max-Cut (plotted 
in figure \12\i 



N 


Median gap 


Error 


16 


0.3203 


0.0056 


24 


0.2323 


0.0057 


32 


0.1844 


0.0057 


64 


0.1113 


0.0058 


128 


0.0496 


0.00473 


160 


0.0291 


0.0049 



There are other possibilities for the scaling of the min- 
imum gap with size, in addition to polynomial or expo- 
nential. For example, we considered a "stretched expo- 
nential" scaling of the form Ae~ cN . Omitting the first 
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Figure 11: (Color online) The gap to the first (even) excited 
state as a function of the adiabatic parameter s for one of 
the N — 128 instances of the Max-Cut problem, showing 
two distinct minima. The first, higher, minimum is close to 
s « 0.36 (the location of the order-disorder phase transition) 
while the other, lower minimum (global in the range) is well 
within the spin-glass phase. 

two points (as we did for the exponential fit) we find the 
fit is satisfactory, as shown in the upper panel of Fig. [13] 
(Q = 0.31). Hence it is possible that the minimum gap 
decreases as a stretched exponential. 

However, if for the instances with more than one mini- 
mum, we just take the minimum close to the critical value 
s c , a different picture emerges, as shown in Fig. 1141 In 
this case, a straight line fit works well for the log-log plot 
(Q = 0.96), but poorly for the log-lin plot (Q = 0.016). 
For consistency, we again omitted the two smallest sizes 
for the log-lin plot. These results indicate that the gap 
only decreases polynomially with size near the quantum 
critical point. 

So far we have plotted results for the median minimum 
gap, which is a measure of the typical value. However, it 
is important to note that there are large fluctuations in 
the value of the minimum gap between instances. This 
is illustrated in Fig. [15] which presents the values of the 
minimum gap for all 47 instances for N = 160. For the 
19 instances with two minima in the range of s studied, 
the minimum at larger s is lower than the one at smaller 
s. For these instances, the figure shows both the "local" 
(smaller-s), and the "global" (larger-s) minima. 

From Fig. [15] we note that a substantial fraction of 
instances for N = 160 have a minimum gap which is much 
smaller than the median, 0.0291(49). Smaller sizes do not 
have such a pronounced tail in the distribution for small 
gaps. We should mention that the gap is not precisely 
determined if it is extremely small because we require 
the condition fiAE > 1. For N = 160 we took j3 = 
2048 so this condition is well satisfied for gaps around 
the median. Hence we are confident that the median is 
accurately determined. However, it is not well satisfied 
for the smallest gaps in Fig. [T31 Thus, while it is clear 
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Figure 12: (Color online) Median minimum gap, on log-linear 
(main panel) and log-log (inset) scales, for the 3-regular Max- 
Cut problem for s in the range 0.3 to 0.5. The straight-line 
fit on the log-linear scale (omitting the two smallest sizes) is 
a much better fit (Q = 0.57) than that of the log-log scale 
(Q = 2.7 x 10 -3 ), in which the two smallest sizes are also 
omitted. 



Figure 14: (Color online) Median minimum gap, on log-linear 
(main panel) and log-log (inset) scales, as a function of prob- 
lem size for the 3-regular Max-Cut problem. Here, the mini- 
mum gaps were taken from the vicinity of the quantum phase 
transition at s = s c ~ 0.36 and are therefore not necessarily 
the global minima. The two fits indicate that in this case the 
polynomial dependence is more probable. 
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Figure 13: (Color online) A stretched exponential fit to 

Ae~ cN for data for the median minimum gap for the 3- 
regular Max-Cut problem, omitting the two smallest sizes. 
The fit is satisfactory (Q — 0.31). 

that a non-negligible fraction of instances for N = 160 do 
have a very small minimum gap, the precise value of the 
very small gaps in Fig. [15] is uncertain. We note that if 
the fraction of instances with a very small minimum gap 
continues to increase with N, then, asymptotically, the 
median would decrease faster than that shown by the fit 



in the main part of Fig. [T^J 



V. SUMMARY AND CONCLUSIONS 

It was demonstrated in Ref. [l9| that the quantum 
adiabatic algorithm fails to solve random instances of 
3-regular 3-XORSAT in polynomial time, due to an ex- 
ponentially small gap in the interpolating Hamiltonian 
which occurs near s = s c = |. This exponentially small 
gap is associated with a first order quantum phase tran- 
sition in the ground state. In this work we have pro- 
vided additional numerical evidence for this. We have 
also demonstrated using a duality transformation that 
the critical value of the parameter s is in fact at exactly 
s c = | . We have also shown that the ground state energy 
of the three regular 3-XORSAT model with a transverse 
field agrees very well with a replica symmetric (RS) cav- 
ity calculation. This provides support for the claim of 
Ref. [l~9| that the RS calculation is exact for the thermo- 
dynamic properties of this model. 

For the random ensemble of Max-Cut instances that 
we consider, we find that the interpolating Hamiltonians 
exhibit a second order, continuous phase transition at a 
critical value s = s c ~ 0.36. Near this critical value of s 
we find that the eigenvalue gaps decrease only polynomi- 
ally with the number of bits. However, we also observe 
very small gaps at values s > s c , i.e. in the spin glass 
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Figure 15: (Color online) A scatter plot of the minimum gap 
for all 47 instances for size N = 160 for the 3-regular Max- 
Cut problem. For the 19 instances with two minima in the 
range of s studied, both are shown, that closer to s = 0.36 
being denoted "local", and the smaller one, at larger s, being 
denoted "global". Note the large scatter in the values of the 
minimum gap for different instances. 

phase. An analysis of the fits indicates that a gap de- 
creasing exponentially with size is preferred over a poly- 
nomially varying gap, though a stretched exponential fit 
is also satisfactory. 

For both of the problems we studied, the adiabatic in- 
terpolating Hamiltonians are stoquastic. This makes it 
possible for us to numerically investigate the performance 
of the QAA using quantum Monte Carlo simulation and 
the quantum cavity method. However it is possible that 
quantum adiabatic algorithms with stoquastic interpo- 
lating Hamiltonians are strictly less powerful than more 
general quantum adiabatic algorithms. 

The QMC calculations consider only instances in which 
the problem Hamiltonian, Hp, has a doubly degenerate 
ground state and a specified value of the ground state en- 
ergy. These instances are exponentially rare. By contrast 
the cavity approach considers random instances. How- 
ever, we do not think that these restrictions invalidate 
the conclusions on the minimum gap summarized in the 
previous paragraph. 

Also, it should be noted that inside the spin-glass 



phase, QMC techniques become less and less efficient as 
the adiabatic parameter s approaches 1, i.e. when the 
Hamiltonian approaches the classical problem Hamilto- 
nian. Hence we have not been able to study s values 
much larger than 0.5 for a broad range of sizes. It is pos- 
sible, indeed likely, that there are other avoided crossings 
in this range which might lead to even smaller minima 
than those found in the studied range, 0.3 < s < 0.5. 
These minima will however not alter the conclusion that 
the overall scaling of the running time of the QAA - when 
applied to the Max-Cut problem - appears to grow expo- 
nentially (or perhaps in a stretched exponential manner) 
with problem size. 

The first order phase transition in 3-regular 3- 
XORSAT prevents the quantum adiabatic algorithm 
from successfully finding a satisfying assignment. In con- 
trast, the second order phase transition in 3-regular Max- 
Cut does not determine the performance of the quantum 
adiabatic algorithm on this problem. In this example 
the small gaps which occur beyond s c cause the quan- 
tum adiabatic algorithm to fail. These small gaps may 
be associated with "perturbative crossings" as described 
in Refs. [IM1- 
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Appendix A: Details of the Quantum Monte Carlo 
Simulations 

In the Max-Cut problem, the global updates are 
achieved by dividing the configurations of the system 
produced within the SSE scheme into clusters and then 
flippin g a fraction of them within each sweep of the simu- 
lation [45| . An important bonus of these cluster updates 
is the existence of "improved estimators" with which one 
considers all possible combinations of flipped and un- 
flipped clusters. Improved estimators are very beneficial 
for determining time-dependent correlation functions as 
the signal to noise is much better than with conventional 
measurements. It is important to note, however, that 
partitioning the SSE configuration into clusters tends 
to be very inefficient as the adiabatic parameter s ap- 
proaches 1, where the entire configuration tends to form 
one big cluster. 

A difficulty arises in extracting the gap for the Max- 
Cut problem due to the bit- flip symmetry of the Hamilto- 
nian for the following reason. Eigenstates of the Hamil- 
tonian are either even or odd under this symmetry (in 
particular, the ground state is even). In the s — > 1 
limit, states occur in even-odd pairs with an exponen- 
tially small gap (see Fig. 2 of Ref. [2l| for an illustration). 
Therefore, the quantity of interest is the gap to the first 
even excited state. We consider correlation functions of 
even quantities, so there are only matrix elements be- 
tween states of the same parity. However, the lowest 
odd level becomes very close to the ground state near 
where the gap to the first even excited state has a min- 
imum. Hence this lowest odd state becomes thermally 
^ populated, with the result that odd-odd gaps are present 
- 1 in the data as well. 



We eliminate these undesired contributions by project- 
4CKg9fflloi]?b^cep3i3Se4Bi45SH£«pa7^®if:^ii%«¥it onian. A 
way of implementing this projection at zero temperature 
is as follows. In standard quantum Monte Carlo simula- 
tions one imposes periodic boundary conditions in imagi- 
nary time t at r = and f3. To project out the symmetric 
subspace one imposes, instead, free boundary conditions 
at r = and /3. The properties of the symmetric sub- 
space can then be obtained, for (3 — > oo, by measurements 
far from the boundaries. We have incorporated this idea 
into the SSE scheme, and use this modified algorithm in 
the simulations of the Max-Cut problem. For the con- 
venience of the reader, this idea is explained in greater 
detail in Appendix |Bj 
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For 3-XORSAT there is no need to employ a projec- 
tion method because the Hamiltonian does not have bit- 
flip symmetry. However, the presence of 3-spin inter- 
actions in the problem Hamiltonian leads to a different 
difficulty. In the Max-Cut case, the (two-spin) Ising in- 
teractions allow the SSE configurations to be partitioned 
into mutually-exclusive clusters which in turn enable the 
use of improved estimators. The nature of the interac- 
tions in the 3-XORSAT problem does not allow for such 
convenient scheme. Rather, construction of clusters for 3- 
XORSAT is done by repeatedly attempting to construct 
clusters using randomly-chosen pairs of spins taken from 
the three-spin operators. The resulting clusters therefore 
have a random component and can in general overlap 
one another. Moreover, it is important to note that not 
all attempts to construct clusters will be successful as in 
some cases the third spin of the operators involved may 
interfere in the construction. In such cases the cluster 
construction must be aborted and restarted. 



for other systems that respect other symmetries if one 
chooses the state |</>) appropriately. 

Next we define a fictitious 'partition function' for the 
scheme, 



Z = (0|0) = (^|e-^/ 2 xe-^/ 2 | 



EE 



Zl ,Z 2 111 ,«2 



1 (P 



ni\n 2 \ \2 



nl+n2 



(B3) 

( Zl \(-Hr(-Hr\ Z2 ). 



As will be immediately clear, the above 'partition func- 
tion' merely serves here as a normalizing factor for the 
various measured quantities. 

As with the usual SSE approach, we divide the Hamil- 
tonian into components, commonly referred to as 'bond' 
operators, 



(B4) 



Appendix B: An SSE-based projection-QMC 
method 



Here we derive in some detail an SSE-based projec- 
tion QMC method. This appendix is rather technical in 
nature and so is intended mainly for readers who are al- 
ready familiar with the SSE method. The purpose of the 
algorithm outlined here is to obtain the zero-temperature 
properties of a system described by a Hamiltonian H pro- 
jected onto an invariant subspace of a discrete symmetry 
operation respected by the Hamiltonian. In this paper, 
we apply the method to N spin-1/2 particles and the 
Max-Cut Hamiltonian, where the goal is to sample only 
states which are even under bit-flip symmetry. It should 
be noted, however, that the method is easily generaliz- 
able to other cases. 

The success of the projection method follows from the 
following observation: Consider the state 



,iV\ 



{*} 



(Bl) 



which is a superposition of all 2 N basis states {z} with 
equal amplitude. Here, the orthonormal set {z} denotes 
the basis of all classical spin configurations along the z- 
direction. Clearly, the state |</>) is symmetric (even) un- 
der bit-flip symmetry. Acting /3 times with the operator 
e~ H / 2 on the state \<f>), where (3 is a large integer, the 
resulting state is, up to a normalization constant, the 
(even) ground state, i.e. 



|0) oc |0) ~e-^/ 2 | 



(B2) 



where |0) is the unnormalized ground state. Note that 
the bit-flip symmetry shared by both the Hamiltonian H 
and the state |</>) confines the projected states to the even 
subspace. This method can be very easily generalized 



The bond operators should obey Hb\z) = h(b, z)\z') for 
all states in {z}. Here, h(b, z) is a real number that 
depends in general on the bond index b and the state z. 
The resulting state \z') must also be one of the 2 N basis 
states chosen for this problem. The partition function 
then becomes 



^=EE E 

»i,»2 ni,n 2 {§ ni ,S n2 } 



1 (fi 



nl+n2 



(zi\S ni S n2 \z 2 ) , 

(B5) 



where S ni (<5 n2 ) denote products or 'sequences' of bond 
operators of length n\ (712) and the summation {S ni } 
({S n2 }) is taken over all possible such sequences. 

We next define a new index variable n = n\ + n%, in 
terms of which the partition function may be written as 

on n 

^EEESE^mn^)- ( B6 ) 



ni=0 



where S n ni ^ denotes an operator sequence of length n 
with an imaginary 'cut' running through it, separating 
the first n\ operators in the sequence from the last (n — 
ni) operators. Here, we have also defined the weights 



w(n,m) = 2~ n 



"1 



which obey 



w(n, n\) = 1 , Vn . 



(B7) 



(B8) 



711=0 



Before moving on, let us first denote by |a(ni)) the (nor- 
malized) state obtained by acting with the first m op- 
erators in the sequence S n on state \z\). In particular 
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|ck(0)) = \zi) and |a(n)) = \z2) ~ Snl^i), where n is the 
number of operators in the sequence. 

The above expression for the partition function 
Eq. (|B6[) has a form very similar to the one obtained 
in the usual SSE decomposition. There are however a 
couple of notable exceptions: (i) While in the usual SSE 
the boundaries in the imaginary time direction are peri- 
odic (the requirement \z\) = \z 2 ) is enforced), here the 
boundary conditions are free and the states \z\) and \z2) 
are different in general and are summed over indepen- 
dently, (ii) To each level along the operator sequence 
there is an assigned weight, reflecting the fact that the 
different time slices in the 'level' direction are not equally 
weighed. The time slices in the middle are weighted more 
than those close to the boundaries. 



1. The updating mechanism 

As in the usual SSE routine, a configuration is de- 
scribed by the pair {|zi), S n }, i.e., a basis state and an 
operator sequence. Importance sampling of the config- 
urations can be done here in much the same way as in 
the usual SSE algorithm. One can use the same local 'di- 
agonal' updating steps by sweeping serially through the 
sequence S n replacing identity operators with diagonal 
ones with appropriate acceptance probabilities and vice 
versa. The acceptance ratios here are exactly the same 
as those in the usual SSE procedure. 

The global non-diagonal updates (normally loop or 
cluster constructions) will also be the same albeit with 
one exception. Here, since the boundaries in the 
imaginary-time direction are free rather than periodic, 
loops or clusters cannot cross the initial and final levels 
to the other side but must terminate at the boundaries. 



2. Static measurements 

The expectation value of an operator A is given by 

(A) = ^ = l(4>\e-^^Ae-^^) (B9) 
= ~z E E ~T E w{n,n 1 ){z l \S ni AS n ^ ni \z 2 ) , 

Z 1' Z 2 71,{S„} ' " 1=0 

where S ni AS n - ril stands for the operator A sandwiched 
between two parts of the sequence S n splitting it in two 
at the cut ri\. The subscripts denote the sizes of each of 
the two sequences, n\ and n — rii, respectively. 

For a diagonal operator, whether a bond operator or 
not, we get 



(i) = 



(BIO) 



1 

Z 



E E £r E Mn,n 1 )a(a ni )](z 1 \S^\z 2 ), 

Z L Z 2 n,{S n } ' ™1=° 

where a(a ni ) — (a(ni)\A\a(ni)) . This means that the 
expectation value of A will be determined from 

n 

(i) = <E Hn,ni)a(a ni )}) . (Bll) 



As in the usual SSE scheme we can only calculate ex- 
pectation values of off-diagonal operators if they are bond 
operators (or products of bond operators) . A general ex- 
pression for the average of a bond operator A (either 
diagonal or off-diagonal) is 



a: 

(A) = Z^ 1 x E E ~ T E w ( n > n i) x (^i|5„ 1 S , Il _„ 1+ i|z 2 ) x S A j W 



711=0 



(B12) 



where 5 ? mi) means that if the first operator in S„_„ +1 is anything other than A then the corresponding weight 

is zero (this is completely analogous to the corresponding derivation in the usual SSE scheme, see, e.g., |36j). Making 
the substitution n — > n — 1, we arrive at: 



on-l n ~ 1 

(A) = Z- 1 X ^ E ( n _ l)\ X E W ( n ~ l , n l){ z l\S ni Sn- ni \ z 2) X J^jW 
n,{£ n } ™ 1=0 



Rewriting the above expression gives 



w(n - 1, ni)—5^ £(i) ) (zi\S ni S n - ni \z 2 ) , 



a> = ^xE E ^r x E 

Z l' Z 2n,{S„} "i=° 

which eventually becomes our final expression for the average of a bond operator: 



ni=0 



w(n, ni) ( (n - "-i)^ s(i) 



711=0 



(B13) 



(B14) 



n—i 

^2 [w(n,ni) ( (n-ni)6 A ^(n-n 1+ i)^ ). (B15) 
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For diagonal bond operators one can use either Eq. (|Blip or (|B15|) . For products of bond operators, we similarly get 

(B16) 



m n—m 



<n^>=(l) x<£ 

i=l VP/ ni=0 



(n-ni)! 

W(n,m) rrd , A(„ I+ l..n 1 + m) 

(n — m— my. iu=i A *> a * 



3. Correlation functions 

For correlation-function measurements, let us consider the following expectation value: 

(A 1 (t/2)A 2 (-t/2)) = {e" T / 2 A 1 e-" T A 2 e" T / 2 ) ra {Q\e" T / 2 A 1 e-" r A 2 e" T / 2 \Q) (B17) 
= (0|ii|0)(0|i 2 |0)+ ^(0|i 1 |m)(m|i 2 |0)e-( £ --^ T . 

m—l 

In our case 

(A 1 (t/2)A 2 (-t/2)) = {e kr ' 2 A^-^ A 2 e" T ' 2 ) ~ (0|e-*^- T ^Aie- AT ii a e~ iW ~ T)A |0) , (B18) 
which becomes 



(i 1 (r/2)i 2 (-r/2)) = Z" 1 x £ £ £ £ y^mX*!^ (ii5 m i 2 ) 5„_„Jz 2 ) . (B19) 

For diagonal operators, this can be expressed as 



(i 1 (r/2)i 2 (-r/2)) = (^ ( /"(I - ^) n " m £ W»-™,«iK(«iW»i + m)] ), (B20) 

m=0 \ ' " " \ni=0 / 

with an analogous expression for bond correlation functions: 

(i,(T/2)i 2 (- T /2)) = (B21) 

(E (;)<£>-<• - §>- Pxf [»(» - m - w A ^ (u .™. ' 

m =o v 7 p p VP/ \ m=0 

4. Integrated susceptibilities 

Integrated susceptibilities are given by 

/ AT{A 1 {T/2)A 2 {-T/2)) = {^—Y j Y j ar{n 1 )a 2 {ni + m)w{n-m,n 1 )), (B22) 

,yu m=0ni=0 

for diagonal operators and by 

~B n— 2 n—m— 2 



ft ^ lb lib ±. 

/ dr(ii(r/2)i 2 (-T/2)) = ^ £ t«(n - m - 2, "O^^i+D^ , 



m=0 ni=0 

for bond operators. 



r 



(B23) 



5. Binomial distribution of the level weights and distribution with p = q = 1/2, the weights are sharply 
the large (3 limit peaked around Hi = n/2, i.e., the mid-point of the se- 

quence. More importantly, the standard deviation of the 
Note that since the weights assigned to the levels, distribution is a = {1/2) meaning that in the limit of 
w(n,m) given in Eq. (|B7). correspond to a binomial very large n, most of the weight is sharply peaked around 
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n/2 and there is no need to perform measurements over 
the entire sequence, as most of the weight is concentrated 
in the region within of order y/ri of n/2. 

We should emphasize, however, that this binomial dis- 
tribution is only needed to reproduce the precise average 
in Eq. (|B10[) for a specific value of /?. However, /3 does not 
correspond to a true inverse temperature and the average 
in Eq. (|B10[) does not, in general, correspond to a Boltz- 
mann distribution. Only for the special case of j3 — >• oo, 
in which limit the method projects out the ground state, 
does this technique give a valid thermal average. In the 
case of large /3 we can, in fact, obtain the ground state 
by sampling anywhere far from the boundaries. For ex- 
ample we can obtain (A) by the following generalization 
of Eq. (|BT0l) . 



(B24) 



where A can take any value between and 1 for which 
both e~ x P H and e~^ x '^ H \<j>) project out the ground 
state. Repeating the above analysis the weights are now 
sharply peaked around m = An. Since different values of 
A give the same result, we can average over A, and hence 
obtain ground state properties by omitting the weights 
w(n,ni) and averaging uniformly over a range of levels 
around the middle (in practice we take the middle n/4 
levels) . Averaging in this way over a range of levels pro- 
portional to n (rather than y/n) improves the signal to 
noise. 



Appendix C: The quantum cavity method for 
two-Local transverse field spin Hamiltonians 

In this section we motivate and describe the equations 
which we have solved numerically in our study of random 
3 regular Max-Cut. We first derive the cavity equations 
for a transverse field spin Hamiltonian with two local 
interactions on a finite tree. We then briefly mention 
the procedure that is used to investigate the infinite size 
limit for homogeneous Hamiltonians defined on random 
regular graphs (homogeneity means that the interaction 
is the same on each edge of the graph). 



1. The quantum cavity method on a tree 

We now review the quantum cavity equations in the 
continuous imaginary time formulation [29j . We consider 
transverse field spin Hamiltonians of the form 



N 



(CI) 



where H is diagonal in the Pauli basis and is two local, 
that is 

Hp = yi H*j 



where only acts nontrivially on spins i and j and the 
graph of interactions T is a tree. 

The starting point for the quantum cavity method 
is the path integral expansion of the partition function 
Tr [e - ^] , where /3 is the inverse temperature. This 
leads to an expression of the form 



Tr 



E 

paths P 



P(P), 



(C2) 



where p is a positive function on paths in continuous 
imaginary time. A path P can be specified by a number 
of flips r a sequence of bit strings {zi, Z2, Z3, 2 r +i = 
zi} where Zi+i differs from zi by a single bit flip, and 
an ordered list of times {ti, £2, ■■■,t r } at which transitions 
occur. By normalizing p we get a probability distribution 
p over paths: 

p(P) = -^-X r dt r dt r . 1 ...dt 1 e-Io(m\Ho\P(t))dt_ 
Z(p) 

A path P of N spins can also be specified as a collection 
of N one-spin paths PW for i e {1, N} where P^ is 
specified by r(i) (the number of transitions in the path of 
the ith spin), a single bit b(i) £ {0, 1} which is the value 
taken by the spin at time t = 0, and a list of transition 
times {ij* , #2 , -,%)}■ Then we can also write 



P(P) 



Z{/3) 



N 



.i=l 



- S$ (P(t)\H \P(t))dt 



(C3) 



The quantum cavity equations allow one to determine 
fii-tj(P( % ') , the marginal distribution of the path of spin 
i when the interaction between spins i and j is re- 



moved from H. 
through 



This marginal distribution is defined 



(pffl) _ 1 p (p-) e tf(p(t)\H t] \p(t))dt^ 



where N^j is a normalizing factor. The quantum cavity 
equations are the following closed set of equations for the 
cavity distributions {pi^j}. 



1 / 



x^dt^d4\..dt% 



e n ^ { p(k: 

p(k) . kedi\j 



k<Edi\j 



■ (P(t)\H, k \P(t))dt 



(C4) 

(C5) 
(C6) 



where 2j_>j is a normalizing constant. From the cavity 
distributions it is straightforward to compute expectation 
values of local operators such as the magnetization or the 
energy. 



19 



2. The thermodynamic limit: replica symmetric 
and 1-step replica symmetry breaking cavity 
equations 

The replica symmetric (RS) scheme is exact under the 
assumption that the measure over paths in Eq. (|C2[) is 
characterized by a single pure state, and local correla- 
tions decay very quickly as a function of distance. In 
this case the loops of the random graph are irrelevant. 
For a model defined on a regular graph, and without dis- 



M(^ (0) ) = -(X r dhdt 2 ...dt r ) 

P(1),P(2) 



One can then attempt to solve for a distribution p over 
paths which satisfies this recursion. Note that if there is 
disorder in the Hamiltonian, e.g. in Eq. (fTi)) . then the 
RS cavity method is more complicated and requires the 
introduction of a distribution of cavity distributions (2(1 
[28| . This is not required in our case. 

The 1RSB ansatz is the next level of refinement within 
the cavity method-it is described in Refs. [46[ and (30| at 
the classical level and in the quantum case in Ref. |4l| . 
The 1RSB ansatz is exact under the assumption that in 
the thermodynamic limit the distribution p in Eq. (|C2[) 
for a random regular graph is a weighted convex combi- 
nation of distributions k which have very little overlap 
(their support is on non-overlapping sets of paths) and 
are uncorrelated. In the 1RSB cavity method the Parisi 
parameter m <G [0, 1] is used to assign the "states" k dif- 



fi{P {0) ) = ^{X r dt 1 dt 2 ...dt r ) 

P(1),P(2) 



is also distributed according to Q. 



Since we cannot represent an arbitrary distribu- 
tion Q(p) in a finite amount of computer mem- 
ory, we represent the distribution Q by a number 
Nd of representatives: that is, marginal distributions 
/ii, /i2j I^n d which are each assigned an equal weight 
in the distribution. Furthermore, each cavity distribu- 
tion fj, is stored as a list of Nr representative paths 
PW, PW, p( NR > which are given weights in the dis- 
tribution w^\w^\...,w (Nn: > (with J2 l w(i) = !)■ 



order in the Hamiltonian, such as the Max-Cut problem 
in Eq. (fT3"j). the local environment of each site is identical 
to all others. Then, in the thermodynamic limit all the 
cavity distributions are the same {pi^j = p for all di- 
rected edges Roughly speaking this assumes that 
a random regular graph is modeled by an "infinite tree" 
which is obtained by assuming translation invariance for 
the recursion in Eq. (|C6|) . For a 3 regular antiferromag- 
net, this gives 



(C7) 

I 

ferent weights in the distribution. By choosing m <G [0, 1] 
appropriately one obtains the correct weighting corre- 
sponding to the distribution p. 

In our study of 3-regular Max-Cut we use the 1RSB 
quantum cavity method with m fixed to be 0. This cor- 
responds to the assumption that each of the distributions 
K is weighted evenly in the distribution p. We made 
this choice here because it greatly simplifies the compu- 
tation [ijj and does not affect much the result for the 
ground state energy of this particular model [27l |. 

We therefore do not present the 1RSB in full 
generality-we now discuss the 1RSB case with m = 0. 
To use this method, we solve for a distribution Q(p) over 
marginal distributions p which has the property: If p\ 
and p,2 are drawn independently from Q, then p defined 
by 



(C8) 

I 

3. Details of the Quantum Cavity Numerics for 
Max-Cut 

Our simulation was run on a Sicortex computer cluster 
in an embarrassingly parallel fashion. We ran two inde- 
pendent simulations at each value of A. We have checked 
our results with a second independent implementation of 
the continuous time cavity method. We used population 
sizes No = 200 and iVp = 15000. We found numerically 
that there is a systematic error associated with taking 
Nr to be too small and that this error increases as j3 is 
increased. We believe that Nr = 15000 is large enough 
to make this error small for our simulation at /5 = 4 (see 



Mi 



p(i) 



112 



p( 2 )) e ~ f£ (P(t)\l<ro<ri+<ro<T2]\P(t))dt 
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(3l| for more details). 

For unknown reasons our computer code sometimes 
(primarily at higher values of the transverse field A and 
larger values of Nr) did not output the data file. This 
computer bug did not seem to compromise the results 



when the output was produced (we checked this by com- 
paring with results from the independent implementa- 
tion). We have only reported data for values of A where 
both independent simulations at Nr = 15000 outputted 
data files. 



